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We describe an iterative method to optimize the multi-scale entanglement renormalization ansatz 
(MERA) for the low-energy subspace of local Hamiltonians on a D-dimensional lattice. For trans- 
lation invariant systems the cost of this optimization is logarithmic in the linear system size. Spe- 
cialized algorithms for the treatment of infinite systems are also described. Benchmark simulation 
results are presented for a variety of ID systems, namely Ising, Potts, XX and Ifeisenberg mod- 
els. The potential to compute expected values of local observables, energy gaps and correlators is 
investigated. 
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I. INTRODUCTION 

Entanglement renormalization [T] is a numerical tech- 
nique based on locally reorganizing the Hilbert space of 
a quantum many-body system with the aim to reduce 
the amount of entanglement in its wave function. It was 
introduced to address a major computational obstacle in 
real space renormalization group (RG) methods [2ll3||4], 
responsible for limitations in their performance and range 



of applicability, namely the proliferation of degrees of 
freedom that occurs under successive applications of a 
RG transformation. 

Entanglement renormalization is built around the as- 
sumption that, as a result of the local character of phys- 
ical interactions, some of the relevant degrees of freedom 
in the ground state of a many-body system can be decou- 
pled from the rest by unitarily transforming small regions 
of space. Accordingly, unitary transformations known as 
disentanglers are applied locally to the system in order 
to identify and decouple such degrees of freedom, which 
are then safely removed and therefore do no longer ap- 
pear in any subsequent coarse-grained description. This 
prevents the harmful accumulation of degrees of freedom 
and thus leads to a sustainable real space RG transforma- 
tion, able to explore arbitrarily large ID and 2D lattice 
systems, even at a quantum critical point. It also leads 
to the multi-scale entanglement renormalization ansatz 
(MERA), a variational ansatz for many-body states [5]. 

The MERA, based in turn on a class of quantum 
circuits, is particularly successful at describing ground 
states at quantum criticality [ll[5l[6l[3[Hl[S[l0l[n] or 
with topological order [12l [13] . From the computational 
viewpoint, the key property of the MERA is that it can 
be manipulated efficiently, due to the causal structure 
of the underlying quantum circuit 0. As a result, it 
is possible to efficiently evaluate the expected value of 
local observables, and to efficiently optimize its tensors. 
Thus, well-established simulation techniques for matrix 
product states, such as energy minimization |14j or sim- 
ulation of time evolution |15j . can be readily generalized 
to the MERA [H [TT] [TH] . 

In this paper we describe a simple algorithm (and sev- 
eral variations thereof) to compute an approximation of 
the low energy subspace of a local Hamiltonian with the 
MERA, and present benchmark calculations for ID lat- 
tice systems. 

Our goal is to provide the interested reader with a 
rather self-contained explanation of the algorithm, with 



enough information to implement it. Sect. [IT] and III 
review and elaborate on the theoretical foundations of 
the MERA jl]]!], and establish the notation and nomen- 
clature used in the rest of the paper. Specifically, Sect. 
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[TTl introduces the MERA, both from the perspective of 
quantum circuits and of the renormalization group, and 
describes several realizations in ID and 2D lattices. Then 
Sect. |III| explains how to compute expected values of 
local observables and two-point correlators. Central to 
this discussion is the past causal cone of a small block of 
lattice sites and the ascending and descending superop- 
erators, which can be used to move local observables and 
density matrices up and down the causal cone. 

Sect. |IV| considers how to optimize a single tensor of 
the MERA during an energy minimization. This opti- 
mization involves linearizing a quadratic cost function for 
the (isometric) tensor, and computing its environment. 
In Sect. |V]we describe algorithms to minimize the energy 
of the state/subspace represented by a MERA. A high- 
light of the algorithms is their computational cost. For 
an inhomogeneous lattice with N sites, the cost scales as 
0{N), whereas for translation invariant systems it drops 
to just O(logiV). Other variations of the algorithm allow 
us to address infinite systems, and scale invariant systems 
(e.g. quantum critical systems), at a cost independent of 
N. 

Sect. |VI| presents benchmark calculations for differ- 
ent ID quantum lattice models, namely Ising, 3-level 
Potts, XX and Heisenberg models. We compute ground 
state energies, magnetizations and two-point correlators 
throughout the phase diagram, which includes a second 
order quantum phase transition. We find that, at the 
critical point of an infinite system, the error in the ground 
state energy decays exponentially with the refinement pa- 
rameter Xi whereas the two-point correlators remain ac- 
curate even at distances of millions of lattice sites. We 
then extract critical exponents from the order parameter 
and from two-point correlators. Finally, we also compute 
a MERA that includes the first excited state, from which 
the energy gap can be obtained and seen to vanish as 
1/N at criticality. 

This paper replaces similar notes on MERA algorithms 
presented in Ref. [TH]. For the sake of concreteness, we 
have not included several of the algorithms of Ref. TS] , 
which nevertheless remain valid proposals. We have also 
focused the discussion on a ternary MERA for ID lattices 
(instead of the binary MERA used in all previous refer- 
ences) because it is somewhat computationally advan- 
tageous (e.g. see computation of two-point correlators) 
and also leads to a much more convenient generalization 
in 2D. 



II. THE MERA 

Let C denote a Z?-dimensional lattice made of N sites, 
where each site is described by a Hilbert space V of finite 
dimension d, so that Yc = V®^. The MERA is an ansatz 
to describe certain pure states |^) S Yc of the lattice or, 
more generally, subspaces Yu C Yc- 

There are two useful ways of thinking about the MERA 
that can be used to motivate its specific structure as a 
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FIG. 1: (Colour online) Quantum circuit C corresponding to 
a specific realization of the MERA, namely the binary ID 
MERA of Fig. |2] In this particular example, circuit C is 
made of gates involving two incoming wires and two outgoing 
wires, p — ptn = Pout ~ 2. Some of the unitary gates in this 
circuit have one incoming wire in the fixed state |0) and can 
be replaced with an isometry w of type (1,2). By making this 
replacement, we obtain the isometric circuit of Fig. |2] 



tensor network, and also help understand its properties 
and how the algorithms ultimately work. One way is to 
regard the MERA as a quantum circuit C whose output 
wires correspond to the sites of the lattice C [5] . Alter- 
natively, we can think of the MERA as defining a coarse- 
graining transformation that maps C into a sequence of 
increasingly coarser lattices, thus leading to a renormal- 
ization group transformation [T]. Next we briefly review 
these two complementary interpretations. Then we com- 
pare several MERA schemes and discuss how to exploit 
space symmetries. 

A. Quantum circuit 

As a quantum circuit C, the MERA for a pure state 
1^1/) G Yc is made of N quantum wires, each one de- 
scribed by a Hilbert space V, and unitary gates u that 
transform the unentangled state jO)**^ into j^*) (see Fig. 

0- . 

In a generic case, each unitary gate u in the circuit C 
involves some small number p of wires, 

u:Y'^P ^ V^f , w^M = uw^ = I, (1) 

where / is the identity operator in V*^^. For some gates, 
however, one or several of the input wires are in a fixed 
state |0). In this case we can replace the unitary gate u 
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Binary ID MERA: 
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FIG. 2: (Colour online) Top: Example of a binary ID MERA 
for a lattice £. with A'^ = 16 sites. It contains two types 
of isometric tensors, organized in T = 4 layers. The input 
(output) wires of a tensor are those that enter it from the top 
(leave it from the bottom). The top tensor is of type (1,2) 
and the rank xt of its upper index determines the dimension 
of the subspace Yu Q Yc represented by the MERA. The 
isometries w are of type (1,2) and are used to replace each 
block of two sites with a single effective site. Finally, the 
disentanglers u are of type (2,2) and are used to disentangle 
the blocks of sites before coarse-graining. Bottom: Under the 
renormalization group transformation induced by the binary 
ID MERA, three-site operators are mapped into three-site 
operators. 



with an isometric gate w 

W -.Win ^"Vout, W^W^l^i^, (2) 

where Vm = Y'^^^" is the space of the pin input wires 
that are not in a fixed state |0) and Yout — V'^^'""' is the 
space of the Pout — P output wires. We refer to u; as a 
(PmyPout) gate or tensor. 

Fig. [2] shows an example of a MERA for a ID lat- 
tice C made of = 16 sites. Its tensors are of types 
(1,2) and (2,2). We call the (1,2) tensors isometries w 
and the (2, 2) tensors disentanglers u for reasons that 
will be explained shortly, and refer to Fig. |2]as a binary 
ID MERA, since it becomes a binary tree when we re- 
move the disentanglers. Most of the previous work for 
ID lattices [2 |5J O [3 UHl [T71 [TH] has been done using the 
binary ID MERA. However, there are many other pos- 
sible choices. In this paper, for instance, we will mostly 
use the ternary ID MERA of Fig. [3] where the isome- 
tries w are of type (1,3) and the disentanglers u remain 



Ternary ID MERA: 
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FIG. 3: (Colour online) Top: Example of ternary ID MERA 
(rank xt, T = 3) for a lattice of 18 sites. It differs from 
the binary ID MERA of Fig. |2] in that the isometries are 
of type (1, 3), so that blocks of three sites are replaced with 
one effective site. Bottom: As a result, two-site operators are 
mapped into two-site operators during the coarse-graining. 



of type (2,2). Fig. |4] makes more explicit the meaning 
of Eq. [2] for these tensors. Notice that describing tensors 
and their manipulations by means of diagrams is fully 
equivalent to using equations and often much more clear. 

Eq. [2]encapsulates a distinctive property of the MERA 
as a tensor network: each of its tensors is isometric (no- 
tice that Eq. [T| is a particular case of Eq. |2]) . A second 
key feature of the MERA refers to its causal structure. 
We define the past causal cone of an outgoing wire of 
circuit C as the set of wires and gates that can affect the 
state on that wire. A quantum circuit C leads to a MERA 
only when the causal cone of an outgoing wire involves 
just a constant (that is, independent of N) number of 
wires at any fixed past time (Fig. |5|. We refer to this 
property by saying that the causal cone has a bounded 
'width'. 

The usefulness of the quantum circuit interpretation of 
the MERA will become clear in the next section, in the 
context of computing expected values for local observ- 
ables. There, the two defining properties, namely Eq. [2] 
and the peculiar structure of the causal cones of C, will 
be the key to making such calculations efficient. 



B. Renormalization group transformation 

Let us now review how the MERA defines a coarse- 
graining transformation for lattice systems that leads to 




past causal cone: 



FIG. 4: (Colour online) The tensors which comprise a MERA 
are constrained to be isometric, cf. Eq. [2] The constraints 
for the isometries w and disentanglers u of the ternary MERA 
can be equivalently expressed (i) diagramatically or (ii) with 
equations. In this paper we will mostly use the diagramatic 
notation, which remains simple for complicated tensor net- 
works. 



a real-space renormalization group scheme, known as en- 
tanglement renormalization [1 . 

We start by grouping the tensors in the MERA into 
T K, log iV different layers, where each layer contains a 
row of isometries w and a row of disentanglers u. We 
label these layers with an integer t = 1, 2, • • • T, with 
T = 1 for the lowest layer and with increasing values of 
T as we climb up the tensor network, and denote by IJr 
the isometric transformation implemented by all tensors 
in layer t, see Figs. [2]and[3j Notice that the incoming 
wires of each IJ-t define the Hilbert space of a lattice Lr 
with a number of sites that decreases exponentially 
with T (specifically, as N2~'^ and A^3~'^ for the binary 
and ternary ID MERA). That is, the MERA implicitly 
defines a sequence of lattices 



-Co Ci 



(3) 



where Cq = C is the original lattice, and where we can 
think of lattice Cr as the result of coarse-graining lattice 

Cr-l. 

Specifically, as illustrated in Figs. [2]and (|3j this coarse- 
graining transformation is implemented by the operator 
Ul that maps pure states of the lattice Ct-i into pure 
states of the lattice Cr, 



(4) 



and that proceeds in two steps (Fig. [6 1. Let us partition 
the lattice Cr-i into blocks of neighboring sites. The 
first step consists of applying the disentanglers u on the 
boundaries of the blocks, aiming to reduce the amount 
of short range entanglement in the system. Once (part 
of) the short-range entanglement between neighboring 




FIG. 5: (Colour online) The past causal cone of a group of 
sites in Co = C is the subset of wires and gates that can affect 
the state of those sites. The example shows the causal cone 
of a pair of nearest neighbor sites of Co for the ternary ID 
MERA. Notice that for each lattice C-^, r = 0, 1,2,3,4, the 
causal cone involves at most 2 sites. This can be seen to be 
the case for any pair of contiguous sites of Cq. We refer to 
this property by saying that the causal cones of the MERA 
have bounded width. 



blocks has been removed, the isometries w are used in 
the second step to map each block of sites of lattice £r-i 
into a single effective site of lattice Cr ■ 

By composition, we obtain a sequence of increasingly 
coarse-grained states. 



l*o> 



l*i> 



1*1 



(5) 



for the lattices {Co,Ci, ■ ■ ■ ,Ct}, where |*i-) = J/^l^r-i) 
and \^o) = \^) is the original state. Overall the MERA 
corresponds to the transformation U = U1U2 ■ ■ ■ Ut, 



(6) 



with 1*0) = t^l^'r). 

Regarding the MERA from the perspective of the 
renormalization group is quite instructive. It tells us 
that this ansatz is likely to describe states with a spe- 
cific structure of internal correlations, namely, states in 
which the entanglement is organized in different length 
scales. Let us briefiy explain what we mean by this. 

We say that the state contains entanglement at a 
given length scale A if by applying a unitary operation 
(i.e. a disentangler) on a region R of linear size A, we 
are able to decouple (i.e. disentangle) some of the local 
degrees of freedom, that is, if we are able to convert the 
state l^*) into a product state |*') (g) |0), where |0) is 
the state of the local degrees of freedom that have been 
decoupled and \'^') is the state of the rest of the system. 
[Here we assumed, of course, that the decoupling is not 
possible with a unitary operation that acts on a subregion 
R' of the region R, where the size A' of R' is smaller than 
the size of R, A' < A]. 

What makes the MERA useful is that the entangle- 
ment in most ground states of local Hamiltonians seems 
to decompose into moderate contributions correspond- 
ing to different length scales. We can identify two be- 
haviors, depending on whether the system is in a phase 
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(11) Ternary MERA scheme for ID lattice 
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FIG. 6: (Colour online) Detailed description of the real-space 
renormalization group transformation for ID lattices induced 
by (i) the binary ID MERA and (ii) the ternary ID MERA. 



characterized by symmetry-breaking order or by topo- 
logical order (see [12] and references therein). In sys- 
tems with symmetry-breaking order, ground-state entan- 
glement spans all length scales A smaller than the corre- 
lation length ^ in the system — and, consequently, at a 
quantum critical point, where the correlation length ^ di- 
verges, entanglement is present at all length scales In 
a system with topological order, instead, the ground state 
displays some form of (topological) entanglement affect- 
ing all length scales even when the correlation length van- 
ishes [11113]. 



C. Choose your MERA 

We have introduced the MERA as a tensor network 
originating in a quantum circuit. Its tensors have in- 
coming and outgoing wires/indices according to a well- 
defined direction of time in the circuit. Therefore, a 
MERA can be regarded as a tensor network equipped 
with a (fictitious) time direction and with two proper- 
ties: 
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(ii] 3x3 MERA scheme for 2D square lattice 
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FIG. 7: (Colour online) Detailed description of the real-space 
renormalization group transformation for a 2D square lattice 
induced by two possible realizations of the MERA, generaliz- 
ing the ID schemes of Fig. [6] In the first case the isometries 
map a block of 2 x 2 sites into a single site, which can be seen 
to imply that the natural size of a local operator, equivalently 
the causal width of the scheme, is 3 x 3 sites. In the second 
case the isometries map a block of 3 x 3 sites into a single site 
and the natural size of a local operator is 2 x 2. As a result, 
the computational cost in the second scheme is much smaller 
than in the first scheme. 
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• Its tensors are isometric (Eq. |2|. 

• Past causal cones have bounded width (Fig. [S]). 

From a computational perspective, these are the only 
properties that we need to retain. In particular, there is 
no need to keep the vector space dimension of the quan- 
tum wires (equivalently, of the sites in the coarse-grained 
lattice) constant throughout the tensor network. Accord- 
ingly, we will consider a MERA where the vector space 
dimension of a site of lattice £,-, denoted Xt, may depend 
on the layer r (this dimension could also be different for 
each individual site of layer t, but for simplicity we will 
not consider this case here). Notice that Xo — d corre- 
sponds to the sites of the original lattice £. 

Bond dimension. — Often, however, the sites in 
most layers will have the same vector space dimension 
(except, for instance, the sites of the original lattice C, 
with xq — d, or the single site of the top lattice Ct, 
with dimension xt)- In this case we denote the domi- 
nant dimension simply by x, and we refer to the MERA 
as having bond dimension x- The computational cost of 
the algorithms described in subsequent sections is often 
expressed as a power of the bond dimension x- 

Rank. — We refer to the dimension xt of the space 
Yct (corresponding to the single site of the uppermost 
lattice Ct) as the rank of the MERA. For xt = 1, the 
MERA represents a pure state j^*) S V^. More gener- 
ally, a rank xt MERA encodes a xr-diniensional sub- 
space Yjj C Yc- For instance, given a Hamiltonian H on 
the lattice £, we could use a rank xt MERA to describe 
the ground subspace of H (assuming it had dimension 
Xt); or the ground state of H (if it was not degenerate) 
and its xt — 1 excitations with lowest energy. The iso- 
metric transformation U in Eq. [6] can be used to build a 
projector P = UU\ 



P:Yc 



P^ = P, tr(P)-XT, (7) 



onto the subspace V[/ C V^. 

Given the above definition of the MERA, many dif- 
ferent realizations are possible depending on what kind 
of isometric tensors are used and how they are intercon- 
nected. We have already met two examples for a ID 
lattice, based on a binary and ternary underlying tree. 
Fig. [7] shows two schemes for a 2D square lattice. It is 
natural to ask, given a lattice geometry, what realization 
of the MERA is the most convenient from a computa- 
tional point of view. A definitive answer to this question 
does not seem simple. An important factor, however, 
is given by the fixed-point size of the support of local 
observables under successive RG transformations — which 
corresponds to the width of the past causal cones. 

Support of local observables. — In each MERA 
scheme, under successive coarse-graining transformations 
a local operator eventually becomes supported in a char- 
acteristic number of sites. This is the result of two com- 
peting effects: disentanglers u tend to extend the sup- 
port of the local observable (by adding new sites at its 



boundary), whereas the isometries w tend to reduce it 
(by transforming blocks of sites into single sites). For 
instance, in the binary ID MERA, local observables end 
up supported in three contiguous sites (Fig. [2]), whereas 
in the ternary ID MERA local observables become sup- 
ported in two contiguous sites (Fig. [3|. 

Therefore, an important difference between the binary 
and ternary ID schemes is in the natural support of local 
observables. This can be seen to imply that the cost of a 
computation scales as a larger power of the bond dimen- 
sion X for the binary scheme than for the ternary scheme, 
namely as O(x^) compared to 0(x*)- However, it turns 
out that the binary scheme is more effective at removing 
entanglement, and as a result a smaller x is already suf- 
ficient in order to achieve the same degree of accuracy in 
the computation of, say, a ground state energy. In the 



end, we find that for the ID systems analyzed in Sect. VI 



the two effects compensate and the cost required in both 
schemes in order to achieve the same accuracy is com- 
parable. On the other hand, in the ternary ID MERA, 
two-point correlators between selected sites can be com- 
puted at a cost 0(x*), whereas analogous calculations in 
the binary ID MERA are much more expensive. There- 
fore in any context where the calculation of two-point 
correlators is important, the ternary ID MERA is a bet- 
ter choice. 

The number of possible realizations of the MERA for 
2D lattices is greater than for ID lattices. For a square 
lattice, the two schemes of Fig. |7]are obvious generaliza- 
tions of the above ones for ID lattices. The first scheme, 
proposed in (see also |5]), involves isometries of type 
(1,4) and the natural support of local observables is a 
block of 3 X 3 sites. The second scheme involves isome- 
tries of type (1, 9) and local observables end up supported 
in blocks of 2 x 2 sites. Here, the much narrower causal 
cones of the second scheme leads to a much better scaling 
of the computational cost with x, only 0(x^^) compared 
to 0(x^^) for the first scheme. 

Another remark in relation to possible realizations con- 
cerns the type of tensors we use. So far we have insisted 
in distinguishing between disentanglers u (unitary ten- 
sors of type p ^ p) and isometries w (isometric tensors 
of type 1 —^ p'). We will continue to use this terminology 
throughout this paper, but we emphasize that a more 
general form of isometric tensor, e.g. of type (2,4), that 
both disentangles the system and coarse-grains sites, is 
possible and is actually used in some realizations [S]. 



D. Exploiting symmetries 

Symmetries have a direct impact on the efficiency of 
computations, because they can be used to drastically 
reduce the number of parameters in the MERA. Impor- 
tant examples are given by space symmetries, such as 
translation and scale invariance, see Fig. [8) 

The MERA is made of 0{N) disentanglers and isome- 
tries. In order to describe an inhomogeneous state \^!) e 
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(i) Non translationallyinvariant MERA {ii)Translationallyinvariant MERA 




(Hi) Scale invariant IVIERA 




FIG. 8: (Colour online) Ternary ID MERA in the presence of 
space symmetries, (i) In order to represent an inhomogeneous 
state/subspace, all disent anglers u and isometrics w are dif- 
ferent (denoted by different colouring). Notice that there are 
N/?> disentanglers (isometrics) in the hrst layer, in the 
second, and more generally in layer r, so that the total 

number of tensors is 2NY}°^^ l/T < 2N . Therefore the 
total number of parameters required to specify the MERA 
is proportional to the size A'^ of the lattice C (ii) In order 
to represent a state/subspace that is invariant under transla- 
tions, we choose all disentanglers and isometrics on a given 
layer of the MERA to be the same. In this case the MERA 
is completely specified by 0(logA'^) disentanglers and isome- 
trics, (iii) In a scale invariant MERA, the same disentangler 
and isometry is in addition used in all layers. 



Yc or subspace V[/ C , all these tensors are chosen to 
be different. Therefore, for fixed x the number of param- 
eters in the MERA scales linearly in N . 

However, in the presence of translation invariance, one 
can use a translation invariant MERA, where we choose 
all the disentanglers u and isometrics w of any given layer 
T to be the same, thus reducing the number of param- 
eters to O(logiV) (if there are T k, \ogN layers). We 
emphasize that a translation invariant MERA, as just 
defined, does not necessarily represent a translation in- 
variant state \^) £ Yc or subspace V[/ C V^. The rea- 
son is that different sites of C are placed in inequivalent 
positions with respect to the MERA. As a result, often 
the MERA can only approximately reproduce translation 
invariant states/subspaces, although the departure from 
translation invariance is seen to typically decrease fast 
with increasing x- In order to further mitigate inhomo- 
geneities, we often consider an average of local observ- 
ables/reduced density matrices over all possible sites, as 
will be discussed in the next section. 

In systems that are invariant under changes of scale, 
we will use a scale invariant MERA, where all the dis- 
entanglers and isometrics can be chosen to be the same 
and we only need to store a constant number of param- 
eters. The scale invariant MERA is useful to represent 
the ground state of some quantum critical systems p] and 
the ground subspace of systems with topological order at 
the infrared hmit of the RG flow [IglfTB] . 



A reduction in parameters (as a function of x) is also 
possible in the presence of internal symmetries, such as 
U{1) (e.g. particle conservation) or SU{2) (e.g. spin 
isotropy). We defer their analysis to Ref. [21]. 

For the sake of concreteness, the explanations in the 
rest of this manuscript refer to the ternary ID scheme of 
Fig. |3] However, analogous considerations also apply to 
any other realization of the MERA. 



III. COMPUTATION OF EXPECTED VALUES 
OF LOCAL OBSERVABLES AND 
CORRELATORS 

Let o^^'^^^^ denote a local observable defined on two 
contiguous sites [r, r-|-l] of C In this section we explain 
how to compute the expected value 

(o["'"+il)v„ =tr(o['^'''+ilp). (8) 

Here P is a projector (see Eq. [?]) onto the XT-dimensional 
subspace Y{j C Yc represented by the MERA. For a rank 
XT = 1 MERA, representing a pure state l^*) G Yc, the 
above expression reduces to 

(o['-'''+il)* = (*|o['-''-+il|*). (9) 

Evaluating Eq. [8] is necessary in order to extract physi- 
cally relevant information from the MERA, such as e.g. 
the energy and magnetization in a spin system. In ad- 
dition, the manipulations involved are also required as 
a central part of the optimization algorithms described 
in Sects. IIVI and [V] The results of this section remain 
relevant even in cases where no optimization algorithm 
is required (for instance when an exact expression of the 
MERA is known [IJIIS]). 

As explained below, the expected value of Eq. (|8| can 
be computed in a number of ways: 

• By repeated use of the ascending superoperator 
the local operator o^'"'''"'"^! is mapped onto a coarse- 
grained operator ot on lattice Ct- Eq. ^ can 
then be evaluated as the trace of the coarse-grained 
operator ot, tr(o['''''"'"^lP) = tr(oT). 

• Alternatively, by repeated use of the descending su- 
peroperator V, a two-site reduced density matrix 
plr,r+i] lattice C is obtained. Eq. ([8| can then 
evaluated as tr(o[''''^+ilP) = tr(o['^^'^+i]p'^''"^^')- 

• More generally, the ascending and descending su- 
peroperators A and V can be used to compute an 

operator Ot and density matrix pr for 

the coarse-grained lattice C^. Eq. (jsj) can then be 

evaluated as tr(o['-''^+iIP) = tr(olr'''' ''''^'')- 

First we introduce the ascending and descending super- 
operators A and V and explain in detail how to perform 



the computation of the expected value of Eq. ([8| . Then 
we address also the computation of the expected value 



{Oh 



tr(OP), 



[r,r+l] 



(10) 



where O is an operator that decomposes as a sum of local 
operators in C; as well as the computation of two-point 
correlators. Finally, we revisit these tasks in the presence 
of translation invariance and scale invariance. 

The ascending and descending superoperators are an 
essential part of the MERA formalism that was intro- 
duced in Ref. [5] (see e.g. Fig. 4 of Ref. 5 for an 
explicit representation of the descending superoperator 
V). These superoperators have been also called MERA 
quantum channel/MERA transfer matrix in Ref. [TU] . 



A. Ascending and descending superoperators 

In the previous section we have seen that the 
MERA defines a sequence of increasingly coarser lattices 
{CqjCi, • • • , J-'t}- Under the coarse-graining transforma- 
tion of Eq. Ul a local operator oJTiTi^^', supported on 
two consecutive sites [r, r-t-l] of lattice £r-i, is mapped 
onto another local operator oV ' supported on two 
consecutive sites [r',r' -1-1] of lattice Ct (Fig. [s]). This is 
SO because in Ulo\.'_-^ Ur most disentanglers and isome- 
trics of Ut and are annihilated in pairs according to 
Eq. [2] The resulting transformation is implemented by 
means of the ascending superoperator A described in Fig. 



El 



(11) 



In order to keep our notation simple, we do not spec- 
ify on which lattice/sites the superoperator A is applied, 
even though A actually depends on r, r and r'. Instead, 
when necessary we will simply indicate which of its three 
structurally different forms (namely, left Al, center Ac 
or right Ar in Fig. [9]) is being used. 

As above, let [r,r -1-1] denote two consecutive sites of 
lattice Cr~i and let [r',r' -1-1] denote two consecutive 
sites of lattice Cr that lay inside the past causal cone of 

[r,r+l] G Cr-i- 

the descending superoperator V of Fig. 



Given a density matrix pi-' in Cr, 



\t t I 1] 

density matrix p\.'_i in £r-i, 



10 produces a 



(12) 



Notice that the descending superoperator T) (which de- 
pends on T, r and r') is the dual of the ascending su- 
peroperator A, V = A* . Indeed, as can be checked in 
Fig. 11 

and Pt , 



by construction we have that, for any o\.Li 



tr(o[Y+^ii?(p[r 



[r y + 1] 



tr 



(13) 




Ascending Superoperators: 



k 



1 w 



± = 




1 r 



zk = 




FIG. 9: (Colour online) The ascending superoperator A trans- 
forms a local operator o^-i of lattice Ct-i into a local opera- 
tor Or of lattice Cr (for simplicity we omit the label [r, r 
that specifies the sites on which Ot~i and Ot are supported). 
Depending on the relative position between the support of 
Or-i and the closest disentangler, the operator can be lifted 
to lattice Ct in three different ways, indicated in the figure as 
(a), (b) and (c). Correspondingly, there are three structurally 
different forms of the ascending superoperator A, namely left 
Ah, center Ac and right Ar, indicated as (a'), (b') and (c'). 
Notice that the figure completely specifies the tensor network 
representation of the superoperator, which is written in terms 
of the relevant disentanglers and isometries (and their Hermi- 
tian conjugates). An e xpli cit form for the average ascending 
superoperator A of Eq. 34 is obtained by averaging the above 
three tensor networks. 



Correspondingly, there are also three structurally differ- 
ent forms of the descending superoperators, namely left 
Vl, center Vc and right Vr in Fig. [lO] 



B. Evaluation of a two-site operator 

We can now proceed to compute the expected value 
tr(o[''.'-+ilp) of Eq. [sjfrom the MERA. This computation 
corresponds to contracting the tensor network depicted 
in the upper half of Fig. 12 
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(b) 



(c) 



Pr-\ 



Pt-1 



Pt-1 



Descending Superoperator: 

(a-) p,., = VM 




k 1 




FIG. 10: (Colour online) The descending superoperator D 
transforms a local density matrix pr of lattice £r into a lo- 
cal density matrix pr-i of lattice jCt-i- Depending on the 
relative position between the support of pr-i and the closest 
disentangler, the density matrix pr canbe lowered to lattice 
jCt-1 in three different ways, indicated in the figure as (a), 
(b) and (c). Correspondingly, there are three structurally dif- 
ferent forms of the descending superoperator D, namely left 
Dl, center Dc and right Vr, indicated as (a'), (b') and (c'). 
An explicit form for the average descending superoperator T) 
of Eq. |40] is obtained by averaging the above three tensor 
networks. 



In a key first step, the contraction of the tensor network 
for tr(o['''''+^lP) is significantly siniphfied by the fact that, 
by virtue of Eq. [2] each isometric tensor outside the past 
causal cone of sites [r, r -t- 1] e C is annihilated by its 
Hermitian conjugate. As a result, we are left with a new 
tensor network that contains only (two copies of) the 
tensors in the causal cone, as represented in the second 
half of Fig. [12] Because the past causal cones in the 
MERA have a bounded width, this tensor network can 
now be contracted with a computational effort that grows 
with N just as O(logA^). One can proceed in several 
ways: 

Bottom-top. — In the bottom-top approach, we 



Duality: 





FIG. 11: (Colour online) The ascending and descending su- 
peroperators, A and T), are dual to each other, see Eq. |13| 
This becomes evident by inspecting the above figure, where 
the superoperators are explicitly decomposed in terms of dis- 
entanglers and isometries. 



would start by contracting the indices of o''"'''^^! and the 
disentanglers and isometries of the first layer (r = 1) of 
the causal cone; then we would contract the indices of 
disentanglers and isometries of the second layer (r = 2) ; 
and so on (Fig. 14). However, this corresponds to 



repeatedly applying the ascending superoperator A on 
Og'^''^^^' = o['''^+^l. Therefore this is precisely how we pro- 
ceed, obtaining a sequence of increasingly coarse-grained 
operators 



[r,r+l] A [ri,ri+l] A [ra.ra + l] A 



ot (14) 



supported on lattices Cq, Ci, £2, and Ct respec- 
tively. Here, the xt x Xt matrix ot at the top of the 
MERA is obtained according to Fig. [13] and the expected 
value of Eq. [Sj corresponds to its trace. 



tr(o 



[r.r+l] p\ ^ 



P) =tr(oT). 



(15) 



Top-bottom. — In the top-bottom approach, we 
would instead start by contracting the indices of the ten- 
sors in the top layer (r = T) of the causal cone; then 
we would contract the indices of the tensors in the layer 
right below (r = T — 1); and so on (Fig. [l5]). How- 
ever, that corresponds to first computing a density ma- 



trix Pt-1 for the two sites of £t-i according to Fig. 13 
and then repeatedly applying the descending superoper- 
ator v. Therefore this is how we proceed, producing a 
sequence of two-site density matrices 



Pt-1 



■ P2 



ra^rs + l] ^[n^n + l] ^{^''''■+11 (16) 



supported on lattices Lt-i, • • • , £-2, Ci and £o respec- 
tively [H]. The last density matrix pI'-^'^+i] = 
describes the state of the two sites of C on which the 
local operator o'''''""'"^! is supported. Therefore we can 
evaluate the expected value of o[''''"+^l , 



tr(o[''''-+ilp) = tr(o[''''^+iV'''''^^^')- 



(17) 



10 



Expected value of a local observable: 




FIG. 12: (Colour online) (i) Tensor network corresponding to 
the expected value tr(ol'''''+yp) of Eq. [s] The two-site oper- 
ator o''"'''"'"^' is represented by a four-legged rectangle in the 
middle of the tensor network. The shaded region represents 
the past causal cone of sites r,r + l £ £.. (ii) All isometric ten- 
sors that lay outside the past causal cone of sites r,r + 1 £ jC 
annihilate and we are left with a simpler tensor network. 



Top-tensor contractions: 




FIG. 13: (Colour online) (i) The top tensor transforms a two- 
site operator ot-i defined on lattice £t-i into a one-site op- 
erator (a XT X Xt matrix) ot on the top of the MERA. (ii) 
The two-site density matrix pr-i on lattice Ct-i is obtained 
through contraction of the top tensor with its conjugate. No- 
tice that pT-i, as well as all pr, are normalized to have trace 
tr(pr) = XT- 



Middle ground. — More generally, one can also eval- 



Bottom-top contraction: 




FIG. 14: (Colour online) The contraction of the tensor net- 
work in the lower half of Fig. [12] using the bottom-top ap- 
proach corresponds to employing the ascending super opera- 
tor A a number of times. In this particular case, we first use 
(i) Ar, then (ii) Ac and then (iii) ^l, to bring the tensor 
network into a simple form whose contraction gives a complex 
number: the expected value of Eq. |8] 



uate the expected value of Eq. |8] through a mixed strat- 
egy where the ascending and descending superoperators 

[r T I 1] 

are used to compute the operator Or ' ^ and density 
matrix Pr ' ^ supported on lattice Cr , which fulfill 

In all the cases above, one needs to use the ascend- 
ing/descending superoperators about T « log N times, 
at a cost O(x^), so that the total computational cost is 
0(x«log7V). 



C. Evaluation of a sum of two-site operators 

In order to compute the expected value 

(0)v, ^tr(OP), 0^^o['-''-+il (19) 
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Top-bottom contraction: 
(i) , (ii) 




FIG. 15: (Colour online) The contraction of the tensor net- 
work in the lower half of Fig. [12] using the top-bottom ap- 
proach corresponds to first implementing (i) a top tensor con- 
traction followed by repeated application of the descending 
super operator D. Specifically, here we first use (ii) "Dl, then 
(iii) T>c and then (iv) ©h, in order to compute the appropri- 
ate density matrix p^^'^+^^ for two sites [r, r + 1] £ £. With 
the density matrix ^['"''"+^1 we can finally compute (v) the ex- 
pectation value tr(o['''''+ilp) = tr(o['''''+ilp''''''"^'')- 



of an operator O on C that decomposes as the sum of 
two-site operators, we can write 



tr(OP) = ^tr(o['-^''+ilp) 



(20) 



and individually evaluate each contribution tr{o^^'^^^^P) 
by using e.g. the bottom-top strategy of the previous 
subsection, with a cost ©(x'^^ log N) . However, by prop- 
erly organizing the calculation, the cost of computing 
tr(OP) can be reduced to 0{x^N'). We next describe 
how this is achieved. The strategy is closely related to 
the computation of expected values in the presence of 
translation invariance, as discussed later in this section. 
Again, there are several possible approaches: 



Bottom-top. 

tors 



We consider the sequence of opera- 



^ 



... Oj 



Oo = O, (21) 



where the operator Or is the sum of local operators. 



N/3'' 



r+1 



(22) 



Ot-_i is obtained from O^-i by coarse-graining. Or ~ 
UlOr-iUr- Each local operator Or in Or is the sum 
of three local operators from Or-i (see (a),(b) and (c) 
in Fig. |9]), which are lifted to Lr by the three different 
forms of the ascending superoperator, Al, Ac and Ar. 
Since Or-i has N/3'^~^ local operators. Or is obtained 
from Or-i by using the ascending superoperator A only 
N/y^^ times. Then, since X^n^o "^^^ < 2, this means 
that the entire sequence of Eq. [21] requires using A only 
O(A^) times. Once Or is obtained, the expected value of 
O follows from 



tr(OP) = tr(OT). 



(23) 



Top-bottom. — Here we consider instead the se- 
quence of ensembles of density matrices 



Et-1 



Ut- 



r El El 

£j2 ^ ^1 



Eq, 



(24) 



where Er is an ensemble of the two-site density 



matrices pr 

Cr, 



[r,r+l] 



supported on nearest neighbor sites of 



Er 



2] J2,3] 
rr 7 



[N/3 



(25) 



From each density matrix in the ensemble Er we can 
generate three density matrices in the ensemble Er-i by 
applying the three different forms of the descending su- 
peroperator, T>L, Vc and Vr (see (a),(b) and (c) in Fig. 
10 1. All the N/y^"^ density matrices of ensemble Er-i 
can be obtained from density matrices of Er in this way. 
Since X]r=o ^ ^' ^^^^ using the descending 

superoperator T) only O(A^) times, we are able to com- 
pute all the density matrices in the sequence of ensembles 
of Eq. [24] Once the ensemble Eq = E has been obtained. 



E ■ 



the expected value of O follows from 
tr(OP) = ^tr( 



(26) 



(27) 



Middle ground. — More generally, we could build 

operator Or as well as ensemble Er and evaluate the 
expected value of O from the equality 



tr(OP) = ^tr(o[r''-+iVl''"+'')- 



(28) 



Each of the strategies above require the use of the 
ascending/descending superoperators 0{N) times and 
therefore can indeed be accomplished with cost O(x^iV). 
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Two-point correlator: 
(i) 



FIG. 16: (Colour online) In order to compute a two-point 
correlator C2{r\,r2) we need to consider the union of the past 
causal cones of sites r\ and r2. Notice that, in contrast with 
the case of a single local operator, the joint causal cone of two 
distant sites typically involves more than two contiguous sites 
of some lattice Cr- This makes the computational cost scale 
as a power of x larger than x*- 



D. Evaluation of two-point correlators 

Let us now consider the computation of a two-point 
correlator of the form 



(29) 



where o^'"! and o^"^ denote two one-site operators applied 
on two arbitrary sites r and s of £, see Fig. |16[ Fig. 
17 shows the tensor network to be contracted. Again, 
we can use Eq. [2] to remove all disentanglers and isome- 
trics that lay outside the joint past causal cone for sites 
r and s. Then, we can proceed to contract the result- 
ing tensor network, for instance through a bottom-top or 
top-bottom approach, with the help of the ascending and 
descending superoperators (and generalizations thereof). 
Notice that since at intermediate layers the two legs of 
the causal cone may contain two sites each one, in gen- 
eral we will need to compute operators/density matrices 
that span more than just two sites, and the cost of their 
computation will be larger than O(x^). 

However, for specific choices of sites r, s e £, we 
are still able to compute C2 {r, s) with overall cost 
O(x^logA^), as illustrated in Fig. 18 We emphasize 



that this was not possible in the binary ID MERA and 
is one of the main reasons to work with the ternary ID 
MERA. For such choices of sites r and s, each of the two 
legs of the joint past causal cone contains just one site 
until, at some layer tq, they fuse into a single two-site 
leg. We can introduce one-site ascendi ng and descend- 
ing superoperators .4*^^-' and 2?^^^ (Fig. 19 1, in terms of 
which we can express, for t < tq, the transformation of 



a product operator oi^^_^ (g) o^^lj^ into a product operator 
oi'-'l ® ol^'^ = A^'Ho^;[,) ® A^'Holt,), (30) 
or of a density matrix p^r ' into a density matrix 

piT:^! = (pw«pa))(p[r^«']), (31) 




FIG. 17: (Colour online) (i) Tensor network to be contracted 
in order to evaluate a two-point correlator C2{ri,r2). Sim- 
ilarly to the case of a local observable Fig. [3] the tensors 
outside of the casual cone annihilate in pairs due to their iso- 
metric character, Eq. [2] The resulting tensor network (ii) is 
much simpler network. However, for a generic pair of sites 
ri,r2 G £,, the joint past causal cone will contain more than 
just two sites per layer, resulting in a computational cost that 
scales with x as a power larger than x^- 



where r, s e Cr-i and r',s' G Cr are sites correspond- 
ing to single-site legs of the causal cone. In, say, the 
bottom-top approach we can compute the correlator of 
Eq. [29] by using the single-site ascending superoperator 
y^'-^' for layers r < tq and then the two-site ascending 
super-operator A for layer t > tq. 



E. Translation invariance 

The computation of the expected value tr(o['''''+^'P) 
of a single local operator o^^'^~^^^ in the case of a trans- 
lation invariant MERA can proceed as explained earlier 
in this section. In the present case one would expect 
the result to be independent of the sites [r, r -I- 1] G C on 
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(<<] 




-c^ 




























• 


1 





an average over sites, 
tr(o['-''-+ilp) 



l^tr(o[-+ilp) 



iV 
1 



(32) 



lj2tr{o^^-^r+np^^-+'^), (33) 



(iii) 




FIG. 18: (Colour online) Two-point correlators for specific 
pairs of sites r, s [at distances of 3' sites for q — 1,2, 3...] can 
be computed with cost 0(x*logA''). This is due to the fact 
that the causal cones for each of r, s contains only one site 
until they meet — (i) at £i, (ii) at £2 or (iii) at £,3. 



(i) o' 



,['■1'] 




Ax,, 



(ii) One-site Ascending 
Superoperator: 

4'i']=^(i)(o[:ii) 



(iii) One-site Descending 
Superoperator: 



Pt-1_i_ 



M = 




FIG. 19: (Colour online) A one-site operator Or-i supported 
on certain sites of jCt-i (corresponding to the central wire of 
an isometry Wr) is mapped onto a single-site operator on £r- 
In this case the (i) ascending and (ii) descending supcropcra- 
tors A^^^ and V^^^ have a very simple form. 



where the terms o^'"'^^'^^ denote translations of the same 
operator o. This average can be computed e.g. by ob- 
taining the N density matrices individually and 
then adding them together, with an overall cost 0{x^N). 
However, with a better organization of the calculation the 
cost can be reduced to 0(x*logiV). 

We first need to introduce average versions of the as- 
cending and descending superoperators. Given a two-site 
operator Or-i in lattice Cr-i, we can build a two-site 
operator Or by using an average of the three two-site op- 
erators resulting from lifting Or-i to lattice Cr, namely 
Al{ot-i), Al{ot-i) and Al{ot-i)- In terms of the av- 
erage ascending superoperator A, 



A = 1^{Al + Ac + Ab), 
this transformation reads 

Or A{Or-l)- 



(34) 



(35) 



Importantly, if we coarse-grain the translation invariant 
operator 



1 \^ Jr,r+1] 



= N/r, (36) 



where N^-i is the number of sites of £t-i and the terms 

o\.j_^ denote translations of O7— 1 , the resulting operator 
can be written as 



where the terms o^^-i'^^ denote translations of Or and 
where Ot and Or-i are related through Eq. 35 In other 
words, the average ascending superoperator A can also 
be used to characterize the coarse-graining, in the trans- 
lation invariant case, of operators of the form of Eq. |19[ 
Let pr denote the two-site density matrix obtained by 
averaging over all density matrices p^T'''^^! on different 
pairs [r, r + 1] of two contiguous sites of Ct, 



(37) 



Nr 



(38) 



which the operator is supported, but a finite bond dimen- 
sion X typically introduces small space inhomogeneities 
in the reduced density matrix p^^'^+^^ and therefore also 
in tr(o[''^'^+ilp) = tr(o['^^''+iV''^'''"^^')- 

Given a two-site operator o, an expected value that is 
independent of [r, r + 1] can be obtained by computing 



and similarly for lattice Cr~i, 



1 ■^-^ [r.r+ll 



Recall that each density matrix on lattice Cr gives 

rise to three density matrices in Cr-i according to the 



(39) 
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three versions of the descending superoperator, namely 
"Dl, T^c and 2?^. It follows that the density matrix pr-i 
can be obtained from the density matrix Pt by using the 

average descending superoperator, 

V=^-{Vl + Vc + Vr), (40) 

that is 

Pr-l=A{pr). (41) 



by using the average ascending/descending superopera- 
tors. This leads to a sequence of operators Or and density 
matrices pr, 

■AAA _ /.„^ 

oo ^ oi ^ • • • ^ ot, oq = o, (48) 

Pa^Pi^ ■■■ ^PT, Pa = P, (49) 

from which the expected value of o is obtained as tr(op), 
as tr(oT) or, more generally, as tr{orPT)- 



We can now proceed to compute the average expected 
value of Eqs. [32]p3l This can be accomplished in several 
alternative ways. 

Bottom-top. — Given a two-site operator o, we com- 
pute a sequence of increasingly coarse-grained operators 



AAA 

Oq ^ Ol ^ 02 ^ 



Ox, 



(42) 



where Or is obtained from Ot-i by means of the average 
ascending superoperator A. Then we simply have 



l^tr(o['--'^+ilp)=tr(oT). 



(43) 



Top-bottom. — Alternatively, we can compute the 
sequence of average density matrices 



Pt-1 



- V _ V _ 
P2 -> Pi -> Po, 



(44) 



where Pr-i is obtained from pr by means of the average 
descending superoperator V and where p = po corre- 
sponds to the average density matrix on lattice £, 



N 



(45) 



in terms of which we can express the average expected 
value as 



^^tr(o['-^'-+ilp) = tr(ap). 



TV 



(46) 



Middle ground. — As costumary, we can also use 
both A and V to compute Or and pr, and evaluate the 
average expected value as 



l5]tr(o['-''-+ilp)=tr(o.p.). 



(47) 



In all the above strategies the average ascending and 
descending superoperators A and V are used 0{log{N)) 
times and therefore the computational cost scales as 
0(x«logA^). 

To summarize, with a translation invariant MERA 
we can coarse-grain a single two-site operator o (with a 
transformation that involves averaging over all possible 
causal cones) or compute the average density matrix p 



F. Scale invariance 

In the case of a translation invariant MERA that is also 
scale invariant, the average ascending superoperators A 
is identical on each layer r, since it is always made of 
the same disentangler u and isometry w. We then refer 
to it as the scaling superoperator S [11 . Its dual S* 
corresponds to the descending superoperator V. 

As derived in Ref. [lOj, the expected value of a local 
observable o in the thermodynamic limit can be obtained 
by analyzing the spectral decomposition of the scaling 
superoperator S, 

5(») = ^ AQ0atr((^Q»), tr(0a0;3) = 5af3- (50) 

a 

We refer to the eigenoperators (j)a of S, 

S{(f>a) = \a(t)a, (51) 

as the scaling operators. Notice that the operators 4>a, 
which are bi-orthonormal to the operators , are eigen- 
operators of 5*, 



S* 



i) — ^a4'a 



(52) 



We recall that the scaling operator S is made of isomet- 
ric tensors (cf. Eq. [2| and therefore the identity operator 
I is an eigenoperator of S with eigenvalue 1 (that is, 5 is 
unital) , 



5(1) = 



(53) 



On the other hand, since the MERA is built as a quantum 
circuit — and descending through the causal cone corre- 
sponds to advancing in the time of a quantum evolution — 
it is obvious that the descending superoperator I? is a 
quantum channel, and so are V and S* (see also [TO])- In 
particular, S* is a contractive superoperator |23j . which 
means that the eigenvalues Aq in Eq. 50 are constrained 
to fulfill |Aq.| < 1. In practical simulations [TP one finds 
that the identity operator I is the only eigenoperator of S 
with eigenvalue one, Ai = 1, and that |Aa| < 1 for a 7^ I. 
Let p denote the corresponding unique fixed point of S* , 



S*{p)^p. 



(54) 



In an infinite system, the local density matrix of any 
lattice Ct (with finite r) results from applying S* on 
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an infinite number of times, and it is therefore equal to 
the fixed point p [10] • Consequently, Eqs. 48 and 49 are 
then replaced with 



5 5 5 
Oo — > Oi 02 03 

, 5* . 5* . 5* . 

p p p p 



Oo = o, 



(55) 
(56) 



where in addition, by decomposing o in terms of the scal- 
ing operators 4>a, 



^Cq(/)q, CQ=tr(0aO), 



(57) 



we can explicitly compute Or'- 

Or = ( 5°-;-o5 )(o) - ^c,,(A„)>„. (58) 



This expression shows that, unless ci ^ 0, the operator 
Or decreases exponentially with r (recall that |Aq| < 1 for 
a ^ I) and its expected value must vanish. The average 
expected value of o then reads: 



lim — Vtr(o[''''^+ilp) = tr(op). 



(59) 



Two-point correlators for selected positions can also 
be expressed in a simple way, by considering the one- 
site scaling superoperator S^^\ wh ich is how we refer to 



the superoperator A^^'' of Fig. 19 



in the case of a scale 



invariant MERA. Its spectral decomposition, 

'S^'^H') = ^Pa1patr{4>a»), tr(Via'0/3) = (60) 

a 

provides us with a new set of (one-site) scaling opera- 
tors V-'q- Given two arbitrary one-site operators o and 
o', we can always decompose them in terms of these ipa 
(similarly as in Eq. 57 1 . Thus we can focus directly on 
a correlator of the form (i/jq'?/'^*'). Here r and s are re- 
stricted to selected positions as in Fig. [18] Then we 
have 



all 



|A„+Afl 



(61) 



where = logg /i^ is the scaling dimension of the scal- 
ing operator tpa, whereas Cap is given by 



Call = tr {{ipa «) V'/3)/0) 



(62) 



with p from Eq. |54[ 

In deriving Eq. [61 [ we have used that, by construction, 
|r — s| = 3'' for some q ~ 1,2,3,---. Coarse-graining 

[rl [si 

ipa ipj3 a number q of times produces a multiplicative 

factor (p.ap.p)'^ and the residual two-site operator ipa'^ip^^\ 
whose expected value gives Cap- On the other hand, by 
noting that ^« = fi^°S!> I'^^^'l = |r - s|1°S3M = |^ _ s|i°S3 
we arrive at (p.aP/i)'^ = |r — s|^°+'^'', which explains the 
denominator in Eq. |61[ 



We note that the polynomial decay of correlators in the 
scale invariant MERA was established in Ref. [5]. Their 
connection with the eigenvalues of the scaling superop- 
erator, Eq. 50 was formalized in Ref. [TD]. Its closed 
expression in Eq. |61[ including the coefficients Ca0, was 
derived in Ref. [llj . where also three-point correlators 
were considered. [Ref. jllj also unveiled a connection 
between the scale invariant MERA and conformal field 
theory. The coefficients Cap of two-point correlators and 
the analogous for three-point correlators are the key to 
identify the operator algebra of primary fields and their 
towers of descendant fields.] 

In conclusion, from the scale invariant MERA we can 
characterize the expected value of local observables and 
two-point correlators, as expressed in Eqs. [59| and [6T][62] 
All critical exponents of the theory can be extracted from 
the scaling dimensions A^. 

The required manipulations include computing p from 
S (using sparse diagonalization techniques) and diago- 
nalizing S^^\ all of which can be accomplished with the 
ternary ID MERA with cost 0(x®). 



IV. OPTIMIZATION OF A 
DISENTANGLER/ISOMETRY 

In preparation for the algorithms to be described in 
the next section, here we explain how to optimize a single 
tensor of the MERA. 

Let H he a, Hamiltonian made of nearest neighbor, 
two-site interactions h^^-'^^^\ 



H 



(63) 



For purposes of the optimization below, we choose each 
term 

hlr.r+l] go i-Jjat it has no positive eigenvalues, 
^[r,r+i] <^ g jxhis Can be achieved with the simple re- 
placement _ Aj^ax-f, where A^ax is the 
largest eigenvalue of ft, ['"'''+11 .] 

Our goal for the time being will be to minimize the 
energy (Fig. 20 i) 



E = tr{HP), 



(64) 



where P is a projector onto the XT-dimensional subspace 
Yjj € Yc , see Eq. ]7] by modifying only one of the tensors 
of the MERA. The optimization of a disentangler u is 
very similar to that of an isometry w, and we can focus 
on describing the latter in more detail. 

Suppose then that, given a MERA, we want to opti- 
mize an isometry w while keeping the rest of the tensors 
fixed. The cost function E is quadratic in w (more specif- 
ically, it depends bi-linearly on w and w^), 



E{w) = tr(^ wMsW^Ns) + ci, 



(65) 



where Mg and are two sets of matrices and Ci is a 
constant (that originates in all the Hamiltonian terms of 
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(i) Energy: f = tr(HP) 




(here Sa>0 are the singular values of T^j)- 

Accordingly, given an initial isometry w, the optimiza- 
tion is performed by iterating the following four steps q„^^ 
times: 

LI. Compute the environment T^, with the newest ver- 
sion of (as explained below, see also Fig 21 1. 



L2. Compute the singular value decomposition 

vswl 



L3. Compute the new isometry w' 
L4. Replace with w'^ . 



-WV^. 



(ii) Isometry Environment: Y,, 




FIG. 20: (Colour online) (i) The energy of a MERA, defined 
E = tv(HP), is represented explicitly in terms of a tensor net- 
work. The removal of an isometry w from this network gives 
(ii) the environment T^, for w (and similarly for disentanglers 
u). By construction we have that E = t^{wTw). 



Eq. 63 outside the future causal cone of w). Unfortu- 
nately there is no known algorithm to solve a quadratic 
problem subject to the additional isometric constraint of 
Eq. [2] One can, however, attempt several approximate 
strategies (see Ref. [TS| for some possibilities). Here we 
describe an iterative approach based on linearizing the 
cost function E{w). 

In this approach, we temporarily regard w and as 
independent tensors, and optimize w while keeping vo^ 
fixed. The cost function reads, up to the irrelevant con- 
stant, simply 



E*{w)=tr{wT^), 



(66) 



where we call the matrix the environment of the isom- 
etry w and we treat it as if it was indepedent of w. E*(w) 
is then minimized by the choice w = —WV^, where V 
and W are the unitary transformations in the singular 
value decomposition of the environment, — VSW^, 



min E*{w) = mmiwVSW^) = -tr(5) = -^8^ (67) 



Enviroments of an isometry: 




FIG. 21: (Colour online) Tensor network corresponding to 
the 6 different contributions to the environment T^, = 
X]i=i '^w of the isometry w. Notice that at each iteration of 
L1-L4 we need to recompute each TJ„ since it depends on the 
updated . Nevertheless, the Hamiltonian term and density 
matrix that appears in remain the same throughout the 
optimization and only need to be computed once. 



The environment of an isometry w (at layer t) 
can be decomposed as the sum of 6 contributions TJ„ 
{i = I,-- - ,6), each one expressed as a tensor net- 
work that involves neighboring isometric tensors of the 
same layer r (disentanglers and isometries) as well as 



one Hamiltonian term h. 



[r,r+l] 



and one density matrix 



\r'r' + l] TT 

Pt , see lig. 
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The two-site Hamiltonian term 



hli^ collects contributions from all the Hamiltonis 
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terms in Eq. 63 included in the future causal cone of the 
sites r, r + 1 of Cr-i and is computed with the help of 
the ascending superoperator A. Similarly, the two-site 

density matrix pr ' Ms computed with the help of the 
descending superoperator P. The computation of /i^l;^ 

and Pt ' , which only needs to be performed once dur- 
ing the optimization of w, has a cost 0(x*log A^). 

On the other hand, once we have h^^^^^^ and pr 
computing T^^ has a cost O(x^) and needs to be repeated 
at each iteration of the steps L1-L4, with a total cost 
O(x^9one)- In actual MERA simulations we find that the 
cost function E^j typically drops very close to the even- 
tual minimum already after a small number of iterations 
(?ono of the order of 10. 

The optimization of a disentangler u is achieved anal- 
ogously, but in this case the environment T„ decomposes 

1,2,3), see Fig. 
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The 



into three contributions T*j (i 
required Hamiltonian terms and density matrices can be 
computed at a cost O(x^logA^), while the optimization 
of u following steps L1-L4 has a cost 0{x^qone)- 



Environments of a disentangler: 




FIG. 22: (Colour online) Tensor networks corresponding to 
the 3 different contributions to the environment T„ = 
X]i=i Til of the disentangler u. Notice that at each iteration of 
L1-L4 we need to recompute each since it depends on the 
updated uK Nevertheless, the Hamiltonian term and density 
matrix that appears in remain the same throughout the 
optimization and only need to be computed once. 



V. OPTIMIZATION OF THE MERA 

In this section we explain a simple algorithm to opti- 
mize the MERA so that it minimizes the energy of a local 
Hamiltonian of the form Eq. |63] We first describe the al- 
gorithm for a generic system, and then discuss a number 



of specialized variations. These are directed to exploit 
translation invariance, scale invariance and to simulate 
systems where there is a finite range of correlations. 



A. The algorithm 

The basic idea of the algorithm is to attempt to mini- 
mize the cost function of Eq. [64] by sequentially optimiz- 
ing individual tensors of the MERA, where each tensor 
is optimized as explained in the previous section. 

By choosing to sweep the MERA in an organized way, 
we are able to update all its 0{N) tensors once with cost 
0{x^N). Here we describe a bottom-top approach where 
the MERA is updated layer by layer, starting with the 
bottom layer r = 1 and progressing upwards all the way 
to the top layer (top-bottom and combined approaches 
are also possible). 

Given a starting MERA and the Hamiltonian of Eq. 
[63) a bottom-top sweep is organized as follows: 



Al. Compute all two-site density matrices p 
all layers t and sites r Cz Ct- 



[r,r+l] 



for 



Starting from the lowest layer and for growing values of 
T=1,2,---,T— 1, repeat the following two steps: 

A2. Update all disentanglers u and isometrics w of layer 



. . . f?' r I 1] 

A3. Compute all two-site Hamiltonian terms hr ' for 
layer t. 

Then, finally, 

A4. Update the top tensor of the MERA. 



In step Al, we compute all nearest neighbor reduced 

density matrices pr ' for all the lattices Cr , so that 
they can be used in step A2. We first compute the density 



matrix for the two sites of Ct-i as explained in Fig. 13 
Then we use the descending superoperator V to compute 
the 6 possible nearest neighbor, two-site density matrices 
of Ct-2- More generally, given all the relevant density 
matrices of Cr , we use V to obtain all the relevant density 
matrices of Cr-i- In this way, the number of operations 
is proportional to the number of computed density ma- 
trices, namely 0{N), and the total cost is 0{x^N). 

Step A2 breaks into a sequence of single-tensor opti- 
mizations that sweeps a given layer r of the MERA. Each 
individual optimization in that layer is performed as ex- 
plained in the previous section. Note that in order to 
optimize, say, an isometry w, we build its environment 
Tm by using (i) the density matrices computed in step 
Al; (ii) the Hamiltonian terms that were either given at 
the start for t = 1 or have been computed in A3 for 
T > 1; (iii) the neighboring disentanglers and isometrics 
within the layer t. We can proceed, for instance, by up- 
dating all disentanglers of the layer from left to right, 
and then update all the isometrics. This can be repeated 
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a number q,^y of times until the cost function does not 
change significantly. 

In step A3, the new disentanglers and isometrics of 
layer r are used to build the ascending superoperator A, 
which we then apply to the Hamiltonian terms of layer 
T — 1 to compute the Hamiltonian terms for layer r. As 
explained after Eq. |22[ each Hamiltonian term in layer r 
is built from three contributions from layer r — 1 . 

In step A4, the optimized top tensor corresponds to 
the XT eigenvectors with smaller energy eigenvalues of 
the Hamiltonian of the two-site lattice Ct-i, obtained 
by exact diagonalization. 

The overall optimization of the MERA consists of it- 
erating steps A1-A4 until some pre-established degree of 
convergence in the energy E is achieved. Suppose this oc- 
curs after q^^ iterations. Then the cost of the optimiza- 
tion scales as 0(x*A^'Zono'7iay'Zitor)- We observe that it is 
often convenient to keep q^^^ and qi^y relatively small (say 
between 1 and 5), since it is not worth spending much ef- 
fort optimizing a single tensor/layer that will have to be 
optimized again later on with a modified cost function. 



B. Translation invariant systems 

When the Hamiltonian H is invariant under transla- 
tions, we can use a translation invariant MERA [2D]. 

In this case, each layer r is characterized by a disentan- 
gler Ut and an isometry Wt- In addition, on each lattice 
Ct we have one two-site hamiltonian /i^. and one aver- 
age density matrix pr- Then a bottom-top sweep of the 
MERA breaks into the steps A1-A4 for the inhomoge- 
neous case above, but with the following simplifications: 

In step Al, we compute Pr-i from _£j. using the average 



descending superoperator V of Eq. 40 



Pr-l = T){pr). 



(68) 



Then the whole sequence {pr-i, ■ ■ ■ ,Pi,Po}, with T w 
log-ZV, is computed with cost 0(x*log A^). 

In step A2, the minimization of the energy E by op- 
timizing, say, the isometry Wt is no longer a quadratic 
problem (since a larger power of Wt appears now in the 
cost function). Nevertheless, we still linearize the cost 
function E and optimize Wr according to the steps Ll- 
L4 of the previous section. Namely, we build the envi- 
ronment (which now contains copies of Wr and wj, 
all of them treated as frozen) , compute its singular value 
decomposition to build the optimal w'^, and then replace 
Wr and wl with w'^. and w'J in the tensor network for T^u 
before starting the next iteration. 

In step A3, the new hamiltonian term h^- is obtained 
from hr-i using the average ascending superoperator A, 



hr-l hr = A{hr-l). 



(69) 



Step A4 proceeds as in the inhomogeneous case. 
The overall cost of optimizing the MERA is in this case 
0(x«log(7V)g„„ 



C. Scale invariant systems 

Given the Hamiltonian H for an infinite lattice at a 
quantum critical point, where we expect the system to 
be invariant under rescaling, we can use a scale invariant 
MERA to represent its ground state [T] [51 [5] [71 [TUl [TT] . 
[A scale invariant MERA is also relevant in the context 
of topological order in the infra-red limit of the RG flow 
[m [13], both for finite and infinite systems; we will not 
consider such systems here]. 

Let us assume that all disentanglers and isometrics are 
copies of a unique pair (u,w). Then, as explained in 
Ref. |llj . the optimization algorithm can be specialized 
to take advantage of scale invariance as follows: 

In step Al, we apply sparse diagonalization techniques 
to compute the fixed point density matrix p from the 
superoperator S* . This amount to applying S* a number 
of times and therefore can be accomplished with cost 
O(x'). 

In step A2, the environment for e.g. the isometry w, 
Til,, is computed as a weighted sum of environments for 
different layers t=1,2,---. In a translation invari- 
ant MERA the environment for layer r is a function 
f{uT,Wr,PT,hr-i) of the pair {ut,Wt), the density ma- 
trix pT and the Hamiltonian term hr-i (specifically, / is 
the sum of the diagrams in Fig. [2T| . A scale invariant 
MERA corresponds to the replacements 



(UrjWr) (U, W), Pt 



P. 



(70) 



so that only /ir-i retains dependence on r. We then 
choose the average environment 



(71) 



where the weight 1 /y reflects the fact that for each isom- 
etry at layer t there are 3 isometrics at layer t — 1 . Using 
linearity of the / in its fourth argument we arrive at 

'^w = f{u,w,p,h), h = ^^hr-i. (72) 

In practice, only a few terms of the expansion of h (say 
r = 1,2,3,4) seem to be necessary. Given T^, the opti- 
mization proceeds as usual with a singular value decom- 
position. 

Steps A3 and A4 are not necessary. 
That is, the algorithm to minimize the expected value 
of H consists simply in iterating the following two steps: 

Sclnvl. Given a pair (u,w), compute a pair (p, h). 

Sclnv2. Given the pairs (w, w) and (p, /i), update the pair 
{u,w). 

In practical simulations it is convenient to include a 
few (say one or two) transitional layers at the bottom 
of the MERA, each one characterized by a different pair 
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Finite range MERA: 




FIG. 23: (Colour online) A finite correlation range MERA 
with T' — 2 layers is used to represent a state of N = 36 sites. 
Since it lacks the uppermost layers, only sites within a finite 
distance or range ~ 3"^ of each other may be correlated. 
More precisely, only pairs of sites (ri,r2) whose past casual 
cones intersect can be correlated, as it is the case with the 
pair of sites (6, 14) but not with (14,26), for which we have 



(ur, Wr)- In this way the bond dimension x of the MERA 
can be made independent of the dimension d of the sites 
of C These transitional layers are optimized using the 
algorithm for translation invariant systems. 



A clear advantage of the finite range MERA is in 
a translation invariant system, where the cost of a 
simulation with range — 3^ is ©(x^loggC), that 
is, independent of N. This allows us to take the 
limit of an infinite system. We find that, given a 
translation invariant Hamiltonian H = ^^j^ /i''^''"^^' , 
where is the same for all r £ C, the op- 

timization of a finite range MERA will lead to the 
same collection of optimal disentanglers and isometries 
{(wi, wi), (u2, W2), • • • , (ut'j wt')}i for different lattice 
sizes N, N' , N" ■ ■ ■ larger than C,. This is due to the ex- 
istence of disconnected causal cones, which imply that 
the cost functions for the optimization are not sensitive 
to the total system size provided it is larger than C,. As 
a result, {(ui, wi), (u2, W2), • • • , (wt', 'U't')} can be used 
to define not just one but a whole collection of states 
\^{N)), \^{N')), \^{N")), for lattices of different 
sizes N,N',N",---, such that they all have the same 
two-site density matrix p and therefore also the same ex- 
pected value of the energy per link. 



{^{N)\h\'i>{N)) = {^{N')\h\^{N')) ^ 



(74) 



D. Finite range of correlations 

A third variation of the basic algorithm consists in set- 
ting the number T' of layers in the MERA to a value 
smaller than its usual one T « log3 N (in such a way 
that the number N^' of sites on the top lattice Ct' may 
still be quite large) and to consider a state j^*) of the lat- 
tice C such that after T' coarse-graining transformations 
it has become a product state. 



l*i> 



1*1 



where |5't') = |0)^^^'. For instance. Fig. 
MERA for TV = 36 and T' = 2. The four top tensors 
of this MERA are of type (0, 3), where the lack of upper 
index indicates that the top lattice Ct' is in a product 
state of its Nt' — 4 sites. 

We refer to this ansatz as the finite range MERA, since 
it is such that correlations in j^*) are restricted to a finite 
range roughly ( Ri 3^ sites, in the sense that regions 
separated by more than ^ sites display no correlations. 
This is due to the fact that the past causal cones of dis- 
tance regions of £ have zero intersection, see Fig. |23] 

Given a ground state j^*) with a finite correlation 
length ^, the finite range MERA with C ~ C turns out 
to be a better option to represent \'^) than the stan- 
dard MERA with T « logg iV layers, in that it offers 
a more compact description and the cost of the simu- 
lations is also lower since there are less tensors to be 
optimized. The algorithm is adapted in a straightfor- 
ward way. The only significant difference is that the top 
isometries, being of type (0, 3), do not require any density 
matrix in their optimization (their environment is only a 
function of neighboring disentanglers and isometries, and 
of Hamiltonian terms) . 
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(73) 



shows a 



In particular, we can use the finite range MERA algo- 
rithm to obtain an upper bond for the ground state en- 
ergy of an infinite system, even though only T' pairs 
{ur,Wr) are optimized. 



VI. BENCHMARK CALCULATIONS FOR ID 
SYSTEMS 



In order to benchmark the algorithms of the previ- 
ous section, we have analyzed zero temperature, low en- 
ergy properties of a number of ID quantum spin sys- 
tems. Specifically, we have considered the Ising model 
[23, the 3-state Potts model [3S], the XX model ^S] and 
the Heisenberg models [27] , with Hamiltonians 



(75) 



r 

iJp„,, = f^^4'^'+ E Mi^Kt-i] (76) 

r \ a=l,2 J 



H 



XX 



-^Heisenberg 



E 



y y 



a Wat 



(77) 
(78) 
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FIG. 24: (Colour online) Top: The energy error of the MERA 
approximations to the ground-state of the infinite Ising model, 
as compared against exact analytic values, is plotted both 
for different transverse magnetic field strengths and differ- 
ent values of the MERA refinement parameter x- The finite 
correlation range algorithm (with at most T = 5 levels) was 
used for non-critical ground states, whilst the scale invari- 
ant MERA was used for simulations at the critical point. It 
is seen that representing the ground-state is most computa- 
tionally demanding at the critical point, although even at 
criticality the MERA approximates the ground-state to be- 
tween 5 digits of accuracy (x ~ 4) and 10 digits of accuracy 
(x = 22). Bottom: Scale-invariant MERA are used to com- 
pute the ground-states of infinite, critical, ID spin chains of 
Eqs. |75||78| for several values of x. In all instances one observes 
a roughly exponential convergence in energy over a wide range 
of values for x ss indicated by trend lines (dashed). Energy 
errors for Ising, XX and Heisenberg models are taken relative 
to the analytic values for ground energy whilst energy errors 
presented for the Potts model are taken relative to the energy 
of a X = 22 simulation. 



where ax , cry and cr^ are the spin 1/2 PauH matrices and 
Mx 1 , Mx 2 and are the matrices 




Mx- 






1 




r 
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, Mx,2 = 1 
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(79) 
(80) 



We assume periodic boundary conditions in all instances 
and use a translation invariant MERA to represent an 
approximation to the ground state and, in some models, 
also the first excited state. For Ising and Potts models the 
parameter A is the strength of the transverse magnetic 
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FIG. 25; (Colour online) The transverse magnetization (az) 
for Ising and | {M^} for Potts, is plotted for translation in- 
variant chains of several sizes A''. Top: For the Ising model, 
the magnetization given from X = 8 MERA matches those 
from exact diagonalization for small system sizes (TV = 6, 18), 
whilst the magnetisation from the A'' = 54 MERA approx- 
imates that from the thermodynamic limit (known analyt- 
ically). Bottom: Equivalent magnetisations for the Potts 
model, here computed with a x = 12 MERA. Simulations 
with larger N systems show little change from the A = 54 
data, again indicating that A = 54 is already close to the 
thermodynamic limit. 



field applied along the z-axis, with Ac = 1 corresponding 
to a quantum phase transition. Both the XX model and 
Heisenberg model are quantum critical as written. 



Fig. 24 shows the accuracy obtained for ground-state 
energies of the above models in the limit of an infinite 
chain, as a function of the refinement parameter x- Sim- 
ulations were performed with either the finite correla- 
tion range algorithm (for the non-critical Ising) or the 
scale invariant algorithm (for critical systems). In all 
cases one observes roughly exponential convergence to 
the exact energy with increasing x- For any fixed value 
of X, the MERA consistently yields more accuracy for 
some models than for others. For the Ising model, the 
cheapest simulation considered (x = 4) produced 5 dig- 
its of accuracy, whilst the most computationally expen- 
sive simulation (x = 22) produced 10 digits of accu- 
racy [28]. The time taken for the MERA to converge, 
running on a 3GHz dual-core desktop PC with 8Gb of 
RAM, is approximately a few minutes/hours/days/ weeks 
for X = 4, 8, 16, 22 respectively. We stress that these sim- 
ulations were performed on single desktop computers; a 
parallel implementation of the code running on a com- 
puter cluster might bring significantly larger values of x 
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0.5 1 
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FIG. 26: (Colour online) Top: Spontaneous magnetization 
{(Tx) computed with a x = 8 MERA for a periodic Ising 
system of TV = 162 sites. The results closely approximate 
the analytic values of magnetization known for the thermo- 
dynamic limit. A fit of the data near the critical point yeilds 
a critical exponent /3mera = 0.1243, with the exact expo- 
nent known as /3ox = 1/8. Bottom: An equivalent phase 
portrait of the Potts model, here with spontaneous magneti- 
zation I (Ma;,i -I- Mj;,2), is computed with a x = 14 MERA 
and is plotted with a fit of the data near the critical point. 
The fit yields a critical exponent /3mera = 0.105 with the 
exact exponent known to be /3cx = 1/9. 



within computational reach. 

Fig. |25] demonstrates the ability of the MERA to re- 
produce finite size effects. It shows the transverse mag- 
netization as a function of the transverse magnetic field 
for several system sizes. The results smoothly interpo- 
late between those for small system sizes and those for 
an infinite chain, and match the available exact solu- 
tions. On the other hand, the MERA can also be used to 
explore the phase diagram of a system. Fig. [26] shows 
the spontaneous magnetization, which is the system's 
order parameter, for a ID chain of iV = 162 sites for 
both Ising and Potts models, where N has been cho- 
sen large enough that the results under consideration do 
not change singnificantly with the system size (thermo- 
dynamic limit) . A fit for the critical exponent of the Ising 
model gives /3mera = 0.1243 whilst the fit for the Potts 
model produces /3mera = 0.105. These values are within 
less than 1% and 6% of the exact exponents /? = 1/8 
and /3 = 1/9 for the Ising and Potts models respectively. 
Obtaining an accurate value for this critical exponent 
through a fit of the data near the critical point is dif- 
ficult due to the steepness of the curve near the critical 
point. Through an alternative method involving the scal- 
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FIG. 27: (Colour online) Two-point correlators for infinite ID 
Ising and Potts chains at criticality (A = 1), as computed with 
X = 22 scale-invariant MERA. Correlators for the Ising model 
are compared against analytic solutions [24] whilst those for 
the Potts model are plotted against the polynomial decay pre- 
dicted from CFT |29]. A scale-invariant MERA produces 
polynomial decay of correlators at all length scales; a fit of 
the form {c^x^ cy'x'^'^^) oc for Ising correlators generated 

by the MERA gives the decay exponent rf" = 0.24996, close 
to the known analytic value 1/4 and similarly for the fits on 
other correlators. Indeed the MERA here reproduces exact 
(cIT'cIT''^'*') correlators for the Ising model at a distance up to 
d = 10^ sites within 0.6% accuracy. Critcal exponents for the 
Potts are also reproduced very accurately. 



ing super-operator S (Sect. Ill I, more accurate critical 
exponents have been obtained in Ref. [11 . 

The previous results refer to local observables. Let us 
now consider correlators. A scale-invariant MERA, use- 
ful for the representation of critical systems, gives poly- 



nomial correlators at all length scales, as shown in Fig. 27 
for the critical Ising and Potts models. Note that Fig. 27 
displays the correlators that are most convenient to com- 
pute (as per Fig.[l8|. These occur at distances d = 3'' for 
(7 = 1,2,3... and are evaluated with cost O(x^). Evalu- 



ation of arbitrary correlators is possible (see Fig. 16 1 but 
its cost is several orders of x more expensive. The pre- 
cision with which correlators are obtained is remarkable. 
A X = 22 MERA for the Ising model gives (crlT'crl''"^'^') 
correlators, at distances up to d = 10^ sites, accurate to 
within 0.6% of exact correlators. Critical exponents 77 are 
obtained through a fit of the form C{r,r + d) oc with 
C as correlators of x, y or z magnetization. For the Ising 
model, the exponents for x,y and z magnetizations are 
obtained with less than 0.02% error each. For the Potts 
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FIG. 28: (Colour online) Top: A x = 8 MERA is used to 
compute the energy gap AE (the energy difference between 
the ground and excited state) of the Ising chains as a 
function of the transverse magnetic field. The gap computed 
with MERA for A'^ = 6, 18 sites is in good agreement with 
that computed through exact diagonalization of the system. 
Inset: Crosses show analytic values of energy-gaps at the crit- 
ical point for — {6, 18,54, 162}. Even for the largest sys- 
tem considered, N = 162, the gap computed with MERA 
Ai^MERA = 9.67 X 10"'^ compares well with the exact value 
AiScx = 9.69 X 10~^. Bottom: Equivalent data for the Potts 
Model where simulations have been performed with a x = 14 
MERA to account for the increased computational difficulty 
of this model. 



model exponents are obtained with less than 0.04% and 
0.08% error for x and z magnetization respectively. 

Finally, we also demonstrate the ability of MERA to 
investigate low-energy excited states by computing the 
energy gap in the Ising and Potts models. Fig. [28| shows 
that the gap grows linearly with the magnetic field A 
and independent of TV in the disordered phase A >> Ac, 
whilst at criticality it closes as 1/N. Even a relatively 
cheap X = 8 MERA reproduces the known critical energy 
gaps to within 0.2% for systems as large as = 162 
sites. The expected value of arbitrary local Hamiltonians 
(besides the energy) can also be easily evaluated for the 
excited state. 



VII. CONCLUSIONS 

After reviewing the conceptual foundations of the 
MERA [TJ [S] (Sect. II-III), in this paper we have pro- 



vided a rather self-contained description of an algorithm 
to explore low energy properties of lattice models (Sect. 
IV-V) , and benchmark calculations addressing ID quan- 
tum spin chains (Sect. VI). 

Many of the features of the MERA algorithm high- 
lighted by the present results can also be observed by 
investigating systems of free fermions [6] and free bosons 
[7] in D — 1,2 dimensions. These include (i) the ability 
to consider arbitrarily large systems, (ii) the ability to 
compute the low-energy subspace of a Hamiltonian, (iii) 
the ability to disentangle non-critical systems completely, 
(iv) the ability to find a scale-invariant representation of 
critical systems and finally (v) the reproduction of ac- 
curate polynomial correlators for critical systems. How- 
ever, the algorithms of Refs. [6l |7] exploit the formalism 
of Gaussian states that is characteristic of free fermions 
and bosons and cannot be easily generalized to interact- 
ing systems. Instead, the algorithms discussed in this 
paper can be used to address arbitrary lattice models 
with local Hamiltonians. 

An alternative method to optimise the MERA is with 
a time evolution algorithm as described in Ref . [T7] . The 
time evolution algorithm has a clear advantage: it can be 
used both to compute the ground state of a local Hamil- 
tonian (by simulating an evolution in imaginary time) 
and to study lattice dynamics (by simulating an evolu- 
tion in real time). We find, however, that the algorithm 
described in the present paper are a better choice when it 
comes to computing ground states. On the one hand, the 
time-evolution algorithm has a time step dt that needs 
to be sequentially reduced in order to diminish the error 
in the Suzuki- Trotter decomposition of the (euclidean) 
time evolution operator. In the present algorithm, con- 
vergence is faster and there is no need to fine tune a 
time step St. In addition, the present algorithm allows 
to compute not only the ground state but also low energy 
excited states. It is unclear how to use the time evolution 
algorithm to achieve the same. 

The benchmark calculations presented in this 
manuscript refer to ID systems. For such systems, 
however, DMRG [3] already offers an extraordinarily 
successful approach. The strength of entanglement 
renormalization and the MERA relies on the fact that 
the present algorithms can also address large 2D lattices, 
as discussed in Ref. [HI W\ ■ 
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